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Abstract -Via simulations of flowing foam, we connect the high and intermediate density regimes 
of complex fluid flows into a consistent microscopic picture of deformation. While at and above 
the jamming transition, elastic correlations lead to strong spatial organization of the flow field, 
below jamming, the slowly diminishing elastic correlation length leads to slowly ceasing spatial 
organization, which is nevertheless still present down to densities far below jamming. We show 
that the long-range correlated flow field arises from the superposition of quadrupolar strain fields of 
shear zones with highly correlated positions, strengths and orientation. These interactions are still 
pertinent below jamming, where they systematically weaken with the slowly diminishing elastic 
correlation length. These results demonstrate the ubiquity and importance of elastic correlations 
in the flow of complex fluids even below the jamming transition, and motivate a scale-bridging 
description of their flow over wide ranges of density from solid to fluid. 


Disordered packings of foams, emulsions, colloidal sus¬ 
pensions and granular particles all exhibit a rigidity tran¬ 
sition to jammed solids when packed densely [iHS]. This 
rigidity vanishes at the jamming transition, where the 
packing approaches its stability limit: the average num¬ 
ber of contacts with neighboring particles approaches a 
critical stability limit, and the shear moduli vanish with 
well-known scaling relations US]. Concomitantly with the 
loss of elasticity, non-afSne fluctuations become increas¬ 
ingly important, and floppy modes indicate the increasing 
susceptibility of the material to applied stress [6HIQ]. At 
unjamming the elastic moduli of the static packing vanish 
and the material is irreversibly affected by the smallest 
applied force [^ISIITO]. 

The situation changes qualitatively when the material 
is subjected to flow. In steady-state flow, particle contacts 
are constantly broken and reformed, leading to a dynamic 
scenario of continuously changing, transient particle con¬ 
tacts m- At the same time, the structural rearrange¬ 
ments have to be relaxed to the boundaries by some long- 
range displacement field. The nature of these fields as the 
density decreases near to and below the static jamming 
transition remains unclear. In the deeply jammed state. 


the material exhibits pronounced elasticity that causes 
correlations in the flow, and provides the long-range field 
that transfers the local relaxation. In this regime, simula¬ 
tions as well as experiments have established that the flow 
of dense foams, emulsions and suspensions is governed by 
local shear transformation zones [HHH] surrounded by 
a long-range quadrupolar elastic strain field [15]. While 
this long-range elastic field - an essential feature of elastic 
materials - provides interactions between transformation 
zones and leads to strongly correlated flow, it is unclear 
what happens near jamming, and how the flow ultimately 
crosses over to the Newtonian regime far below jamming. 
Despite its central importance for many applications of 
complex fluids in industry and consumer products, this 
intermediate flow regime remains poorly understood. 

In this letter we elucidate just this intermediate flow 
regime between the strongly correlated high-density, and 
the low-density Newtonian regime using simulations of 
flowing foam. We show that while similar to static pack¬ 
ings, the flow field becomes increasingly delocalized with 
decreasing density, it is still uniquely determined by elas¬ 
tic correlations at and below the jamming density, where 
the elasticity of the static packing vanishes. The robust 
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Fig. 1: Non-afBne displacements during the flow of foam at 
shear-rate 7 = 10 ~® for bubble volume fractions of 0 = 1.0 (a) 
and (j> ^ (j>j = 0.8424 (b). The flow field becomes delocalised. 


power-law correlations persist across the jamming tran¬ 
sition and become systematically truncated by a slowly 
diminishing cut-off below jamming. We demonstrate that 
the non-trivial power-law exponent arises from the super¬ 
position of correlated quadrupolar strain fields that due 
to ceasing elastic interactions become increasingly delocal¬ 
ized. These results demonstrate the abundance of robust 
elastic-like correlations even in regimes where the elastic¬ 
ity of the static packing vanishes, highlighting the crucial 
role of elastic correlations in complex fluid flows, and mo¬ 
tivating a universal scale-bridging framework over a wide 
range of density. 

Model: We probe correlations of flowing foam in 
Durian’s two-dimensional bubble model m, which cap¬ 
tures many aspects of complex fluid flows above, near, and 
below jamming mnn- In this model, particles or bubbles 
are represented by disks which interact through purely re¬ 
pulsive elastic and viscous contact forces. Inertia is absent, 
and the sum of elastic and dissipative forces on each parti¬ 
cle balance at all times. Elastic forces are proportional to 
the disk overlap/j®* = k{Ri+Rj—rij), where Vij := \fj—fi\ 
is the distance between bubble centers and Ri is the radius 
of disk i. Viscous forces oppose the bubbles’ relative ve¬ 
locity Avij := Vj — Vi with magnitude = b\Avij\. 

The packing fraction </) controls the density, with the 
jamming density of order t/fj = 0.8424 [18]. The strain 
rate 7 is imposed via Lees-Edwards boundary conditions, 
which lead to a linear flow profile where (u(t/)) = 72 / e^,, 
with y and x the transversal and stream-wise coordinates. 
The unit cell contains a 50 : 50 bidisperse mixture of 
N = 1020 — 1210 bubbles with size ratio 1.4 ; 1 to avoid 
crystallization. In the following we use particle diameter, 
elastic (fc) and viscous ( 6 ) prefactors to nondimensionalize 
our results. The global shear stress is obtained according 
to (J^y = atot = l/( 2 E)E<y>nj>(/f/y+/" 7 y°), where V is 
the area of the unit cell and the sum runs over contacting 
pairs. We focus here on slow shear, namely foam sheared 
at a rate of 7 = 10“®. Averages are taken over runs of 
total time 2 O /7 after discarding the transient. 

To probe fluctuations in the local strain field, we fol¬ 
low [Hdi] and determine, for each particle, the local 
strain from the displacement of a particle with respect 





Fig. 2: Strain correlations and scaling. (a,b) Spatial correla¬ 
tions of strain for flowing foam with 0 = 1.0 (a) and 4> = 0.740 

(b) . Color indicates valne of the normalized correlation func¬ 
tion, see color bar. Symmetry change indicates loss of elastic¬ 
ity. (c,d) Decay of strain correlations and correlation length. 

(c) Radial decay of strain correlations projected onto the 4th 
circular harmonic. Correlations are power-law for (j) > (f)j, 
and become increasingly short-ranged below tjjj (d) Correla¬ 
tion length (^) extracted from the decay of correlations in (c), 
as a function of distance A(j> from jamming. Error bars are 
twice the symbol size. 


to its nearest neighbors. We identify nearest neighbors as 
those separated by less than ri, the radius of the larger 
bubbles. We subsequently determine the best affine de¬ 
formation tensor F that transforms the nearest neighbor 
vectors, di, over the applied strain interval [SI, by min¬ 
imizing = (l/n)X;r=i(di(7 + < 57 ) - rdi(7))2. The 

symmetric part of F is the local strain tensor, whose off- 
diagonal component is the shear strain e = which we 
will use in the further analysis below. 

Phenomenology: Snapshots of the deformation field 
with mean flow subtracted reveal significant spatial and 
temporal heterogeneities, and moreover, a distinct trend in 
their qualitative nature with packing density (j), as shown 
in Fig [Tj Well above the jamming density, the deforma¬ 
tion field appears strongly localized with clear quadrupo¬ 
lar structures (Fig.[T^), while at densities approaching <f)j, 
rearrangements become more extended and the displace¬ 
ment field becomes delocalized (Fig. |T|d). Both the strong 
spatial fluctuations and trends with packing density re¬ 
semble those of quasistatic deformations, as studied ex¬ 
tensively in the context of jamming BM- 

Correlations: The question we wish to answer is, what 
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do these fluctuations tell us about the underlying nature of 
the deformation field? In particular, how does the high- 
density physical picture of local rearrangements coupled 
through an elastic background change when we approach 
the jamming point, or even go to densities below jamming? 

To probe the underlying spatial organization of the in¬ 
stantaneous deformation fields, we focus on the local shear 
strain e, and compute spatial correlations according to 
[H] Ce{Sf) = {{e{f + 6 f)e{f)) — (e(f))^)/cr^, where the 
squared standard deviation = (e(f)^) — (e(f))^. As 
shown in Fig. [2^ and b, these correlations reveal the typi¬ 
cal symmetries and range of strain fluctuations. For large 
densities (Fig. [2^), Ce{Sx, Sy) has a distinct four-fold sym¬ 
metry stemming from the abundance of local quadrupolar 
strain fields that reflect the response of an elastic matrix 
to local shear transformations m- Surprisingly and as we 
will show in detail below, this quadrupolar elastic response 
remains essentially unchanged down to the jamming tran¬ 
sition, where the instantaneous deformation field does not 
any more show clear quadrupoles (Fig. lb). In fact, the 
quadrupolar symmetry of C^{Sx, Sy) vanishes only far be¬ 
low jamming, where the correlation function eventually 
becomes isotropic and short ranged (Fig. [^b)- The ob¬ 
servation that elastic correlations behind this flow field 
remain important near jamming is also consistent with ob¬ 
servations of system-spanning velocity correlations of the 
instantaneous flow field [TTIIT^ . 

To quantify the gradual loss of elasticity, we investigate 
the range of quadrupolar correlations. We project each 
correlation function onto the corresponding circular har¬ 
monic, using Cl{Sr) = Ce(Sr,0)cos(4O)d0, and study 
the radial decay of C 4 . In Fig. He) we show examples 
of Cii^Sr) for densities above, near and below jamming. 
These show that down to 4>j, the quadrupolar field is long 
ranged: strain correlations remain robust upon approach¬ 
ing the jamming point. Below jamming, the range of corre¬ 
lations decreases slowly until at ^ ^ 0.74, the correlation 
length has decreased down to a few particle diameters, 
and the correlation function appears closely isotropic. We 
extract a correlation length quantitatively by fitting the 
decay of correlations to Cl 5r°‘ exp(—dr/^) with the 
exponent a = —0.7, yielding excellent overlap with the 
measured correlation function. In Fig. H'i) plot the 
resultant correlation length ^ as a function of density dif¬ 
ference Acj} = (j) — (jij, focussing on At/) < 0, i.e., below 
jamming. We find that this length scale appears to di¬ 
verge as ^ ^ |A(/)|“^, with the best fit for /? = 0.3 ± 0.1. 

Quadrupoles: The powerlaw radial decay of Ci evi¬ 
dences strong organization of the flow in the jammed 
and near jammed regime, irrespective of At/). What are 
the basic ingredients that generate this flow field? The 
quadrupolar correlation function in Fig. [2^ suggests that 
the strain field is composed of an abundance of local 
quadrupolar strain fields, each of which is the elastic re¬ 
sponse to a local shear rearrangement. As we will show, 
the connection between quadrupoles and C 4 is subtle but 




Fig. 3: (a) Synthetic strain field constructed for a single Es- 
lielby inclusion. We selected a particle with large shear strain 
(e = 0.5) as center of the inclusion, and computed the shear 
strain at the positions of all other particles using eq.[T] (b) De¬ 
cay of the projected strain correlations computed for the single 
Eshelby inclusion shown in (a). 



firm. First, we note that a single Eshelby quadrupolar in¬ 
clusion decays as in contrast to C 4 which decays as 

Sr°‘ with a ~ —0.7 > —2. Second, we will show that re¬ 
alistic deformation fields are composed of a superposition 
of many of such quadrupoles, and that a superposition of 
quadropoles, precisely placed at the locations obtained in 
the simulations, reproduces a powerlaw radial decay of C 4 
with an exponent Ri —I > —2. Third, a superposition 
of the same quadrupoles placed at random positions does 
lead to C 4 ~ I/dr^, a similar decay as for the single inclu¬ 
sion. Taken together, this shows that a nontrivial spatial 
organization of quadrupoles is a necessary and also suf¬ 
ficient ingredient to reconstruct the nontrivial powerlaw 
radial decay of C 4 . 

We start from a single quadrupole. Since (weakly) 
jammed materials behave on average as a linear elastic 
matrix [20], we use the elastic shear strain produced by a 
single sheared inclusion of size a, with shear strain eg. 


€^{Sr, 0) 


eg cos(4(6> - dp)) 
TT (dr/a)^ 


Sr » a, 


( 1 ) 


where 0g denotes the orientation of the quadrupole [la 
Unless noted otherwise, we take 0g = 0. Clearly, 
decays as Sr~^, which is qualitatively different from 
the decay of the correlation functions — so a single Es¬ 
helby quadrupole cannot capture the radial decay of C 4 . 
Before studying superpositions of , we show in Eig. 3 
an example of a single quadrupolar field constructed by 
applying Eq. O to a given packing. Strictly speaking, 
Eq. Il]) is only valid at sufficient distance from the center 
of the sheared inclusion, but here we use this form all the 
way to a central particle with large shear strain (e = 0.5), 
and compute the shear strain at the positions of all other 
particles using Eq. O- For the center particle, we take 
the actual shear strain of the particle as determined from 
the affine fitting; repeating this procedure for ^ 1000 re¬ 
alizations (separated by strains of ^ 5 • 10“^) we obtain 
the corresponding correlation function calculated for sin¬ 
gle, isolated shear events, which clearly decays as 1/Sr^ 
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Fig. 4: Correlations of superimposed strain fields, (a) Snapshot 
of the instantaneous strain field and selected inclusion centers 
(dotted particles) at <() = 1.0. Particle color indicates sign 
and magnitude of the local shear strain {txy), see color bar. 
(b) Strain correlations as a function of 5r after quadrupolar 
projection, see text. Data taken from simulations at volume 
fractions 1.0 and 0.8424 (red line), from superimposed inclu¬ 
sions positioned at strain maxima (red circles), from inclusions 
at random positions (blue squares), and from inclusions at the 
right positions, but with random orientations (black asterisks). 
The superimposed inclusions describe the strain field correctly 
if they are at the right positions. (c,d) Crossover upon varia¬ 
tion of number of inclusions and particle volume fraction, (c) 
Dependence of constructed strain correlations on the number 
of inclusions. Taking a few top-strain inclusions leads to the 
trivial decay with slope -2, and does not capture the real cor¬ 
relations. (d) Volume fraction dependence of pair correlations 
of inclusions, see legend. With decreasing volume fraction, 
and concomitant decreasing elastic correlation length, the spa¬ 
tial distribution of inclusions approaches that of the randomly 
placed inclusions. 


(Fig.Eb). 

Synthethic Strain Fields: We now turn to synthetic 
strain fields that we create by superposing many Eshelby 
quadrupoles given by Eq. ©■ We show a snapshot of the 
instantaneous strain field for cj) = 1.0 in Fig. 0^; clearly, 
many locations of significant strain play a role here. We 
thus constructed strain fields for 50 to 70 inclusion cen¬ 
ters identified as particles having shear strain amplitudes 
larger than two times the standard deviation (dotted par¬ 
ticles in Fig.H^). We used eq.[T] where we took eo to be the 
actual strain of the particle at the inclusion center, and a 
to be the particle diameter. We then computed the total 
strain field by superposition, and determined strain cor¬ 


relations and quadrupolar projections as before. The re¬ 
sulting correlation function, shown in Fig. 4b (red circles), 
exhibits excellent overlap with the one obtained in the sim¬ 
ulations (red line). Hence, synthetic strain fields obtained 
by superposition of a sufficient number of quadrupoles can 
capture the nontrivial power law decay of C^. 

Which microscopic properties of the inclusions are es¬ 
sential to capture these nontrivial correlations? We first 
tested what happened when we randomize the locactions 
of the quadrupole centers. We no longer chose high- 
strain particles, but instead selected particles randomly 
(again between 50 and 70) as core of inclusions; we then 
constructed the resulting strain field by superposition, 
and computed the corresponding correlation functions, see 
Fig. 0 )d, blue squares. Clearly, the data does not match 
the nontrivial decay of the correlations and rather decays 
with slope -2, similar to the single inclusion. This shows 
that a nontrivial spatial organization of the inclusions is 
essential to obtain the nontrivial power-law decay of C^. 

We further tested the robustness of these results with re¬ 
spect to the orientation of the inclusions. Typically, inclu¬ 
sions align their directions along the principal directions 
of the system to minimize their interaction energy, as seen 
in Fig. [2^. At finite shear rate, however, we may expect 
some randomness in the orientation of the inclusions. We 
therefore added a random offset angle Oq-, uniformly dis¬ 
tributed between —7r/4 and 7r/4 in Eq. [T] The resulting 
projected strain correlation (asterisks in Fig. 0Kb)) shows 
nice overlap with the power-law data of the real inclusions 
as long as the angle variation is not too large, indicating 
that with respect to the orientations, the scaling of corre¬ 
lations is robust. 

We also investigated variations in the number of inclu¬ 
sions. Clearly, using one inclusion does not work, but a 
sufficient number of inclusions robustly pushes the power- 
law decay away from a slope -2, as illustrated in Fig. 4c. 
Hence, a broad population of inclusions is needed to cap¬ 
ture the full nontrivial behavior. Similarly, when we nar¬ 
row strain amplitudes eo to one fixed value for all inclu¬ 
sions, the resulting correlations weakly deviate from the 
-1 slope, although the effect is not as strong as for ran¬ 
domizing the locations of the quadrupoles. We conclude 
that the inclusions are strongly correlated in their posi¬ 
tion, with less strong correlations in their strength and 
orientation. These correlations are directly reflected in 
the spatial organization of the quadrupoles. To show this, 
we have computed pair correlations of the inclusion cen¬ 
ters [23]; at (/> = 1, we observe much stronger short-range 
correlation than for the random arrangement, as shown in 
Fig.0Jl. 

How does this picture change below (f>c^ Below 
jamming, elastic correlations become increasingly short 
ranged, and this should affect the organization of inclu¬ 
sions. This is indeed what we find when we investigate 
the inclusion pair correlations as a function of volume frac¬ 
tion across jamming, see Fig.0Ji. As the density decreases, 
the pair correlation function becomes increasingly flat, re- 
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fleeting the diminishing spatial organization of the inclu¬ 
sions, until dX 4> = 0.74, the pair correlation function ap¬ 
proaches the randomly placed inclusions (while still being 
distinct from it), suggesting only small remaining inter¬ 
actions, in agreement with the elastic correlation length 
shown in Fig. [51 

Conclusions: In conclusion, we have investigated the 
nature of the flow of foams at hnite but small shear rate 
in a numerical model. When the packing fraction is well 
above jamming ((j) > (j)c), isolated plastic events dominate 
the deformation. These plastic events self-organize due to 
their quadrupolar elastic interactions. Surprisingly, these 
elastic interactions remain robust and long-range down to 
the jamming transition, where the elastic modulus of the 
quiescent packing vanishes, and the displacement field be¬ 
comes delocalized. Thus, down to and even below the jam¬ 
ming density, the spatial organization of flowing foams is 
dominated by long-range elastic interactions. This is con¬ 
firmed by modeling of the strain field correlations using the 
well-known Eshelby elastic response of a single shear dis¬ 
tortion. We showed that the nontrivial flow field leads to 
the superposition of strain helds at strongly correlated lo¬ 
cations. These results shed new light on the nature of flow 
across the jamming transition: over wide ranges of den¬ 
sity from the liquid via the marginal solid to the deeply 
quenched solid state, complex fluid flows arise from the 
superposition of elastic fields of interacting Eshelby inclu¬ 
sions, with interaction ranges only becoming hnite below 
the jamming transition. 
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